1 Setup

1.1 Libs, directory

library("stringr")
library("dplyr")
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library("tidyr")
library("cowplot")
library("ggplot2")
library("cowplot")

setwd("~/nicolagk@hawaii.edu - Google Drive/My Drive/Bamboo_mesos/Bamboo_mesos/04.mixed_models")

1.2 Data

#larv.ranef <- read.csv("larv.mod.ranef.csv")
larv.ranef <- read.csv("larv.mod.ranef.all.csv")

##ps object & taxa table
#ps.lar.trim <- readRDS("ps.exp.lar.trim.rds")
ps.lar.trim <- readRDS("ps.exp.all.lar.trim.rds")
tax.lar.trim <- data.frame(ps.lar.trim@tax_table)
## Loading required package: phyloseq
#tax.lar.trim

1.3 Checking raw stuff

##individual ASVs in all treats
##raw counts
ps.lar.trim.asv <- subset_taxa(ps.lar.trim,id=="ASV0001")
plot_bar(ps.lar.trim.asv)+
  facet_wrap(Food_type~Microbe_treatment,scales="free",ncol=5)

ps.lar.trim.rel = transform_sample_counts(ps.lar.trim, function(x) x / sum(x))
#ps.lar.trim.rel.trt <- subset_samples(ps.lar.trim.rel,Microbe_treatment=="MEM")
ps.lar.trim.rel.asv <- subset_taxa(ps.lar.trim.rel,id=="ASV0001")

plot_bar(ps.lar.trim.rel.asv)+
  facet_wrap(Food_type~Microbe_treatment,ncol=5,scales="free")

1.3.1 Reformatting ranef table

table(larv.ranef$grpvar)
## 
##                             variable                   variable:Food_type 
##                                   54                                  108 
##           variable:Microbe_treatment variable:Microbe_treatment:Food_type 
##                                  270                                  540
#table(larv.ranef$term)

larv.ranef.treat <- subset(larv.ranef,grpvar=="variable:Microbe_treatment:Food_type")

larv.ranef.treat$id <- substr(larv.ranef.treat$grp,1,7)
larv.ranef.treat$food <- str_sub(larv.ranef.treat$grp,13,18)
larv.ranef.treat$microbes <- str_sub(larv.ranef.treat$grp,9,11)
larv.ranef.treat$food_micr <- paste0(larv.ranef.treat$food,"_",larv.ranef.treat$microbes)
#larv.ranef.treat$food_micr <- substr(larv.ranef.treat$grp,9,18)

##just the columns for stuff below:
larv.ranef.less <- larv.ranef.treat %>%
  select(-X,-term,-grp,-condsd,-food,-microbes)

# Reshape the data frame to have the food types as two columns
larv.ranef.wide <- larv.ranef.less %>%
  pivot_wider(
    names_from = food_micr,   # Column to use for new column names
    values_from = condval       # Column to fill the new columns
  )

##add taxa info
larv.ranef.wide.tax <- merge(larv.ranef.wide,tax.lar.trim,by="id")
##more informative ASV ids
larv.ranef.wide.tax$asv.gen <- paste0(larv.ranef.wide.tax$id,"_",larv.ranef.wide.tax$Genus)

2 Uncorrected conditional modes

library("pheatmap")
library("RColorBrewer")

cond.df <- larv.ranef.wide.tax[,c("asv.gen","Bamboo_EMO","Bamboo_MMO","Bamboo_MEM","Larval_EMO","Larval_MMO","Larval_MEM")]

#subset(cond.df)

row.names(cond.df) <- cond.df$asv.gen
cond.df2 <- cond.df[,2:7]
cond.mat <- as.matrix(cond.df2)

#row.names(otu.rel) == row.names(sam.rel)

pheatmap(cond.mat,
         #cluster_cols=F,
        #color = colorRampPalette(c("blue", "white", "red"))(100),
        color = colorRampPalette(rev(brewer.pal(n=7,name="PiYG")))(100)
)

pheatmap(cond.mat,
         scale="row",
         #cluster_cols=F,
        #color = colorRampPalette(c("blue", "white", "red"))(100),
        color = colorRampPalette(rev(brewer.pal(n=7,name="PiYG")))(100)
)

3 Conditional modes by treatments

3.1 Calculate differences

3.1.1 Between foods

##MMO - MMO BETWEEN FOODS
larv.ranef.wide.tax$diff.bmmo.lmmo <- larv.ranef.wide.tax$Bamboo_MMO - larv.ranef.wide.tax$Larval_MMO
##EMO - EMO BETWEEN FOODS
larv.ranef.wide.tax$diff.bemo.lemo <- larv.ranef.wide.tax$Bamboo_EMO - larv.ranef.wide.tax$Larval_EMO
##MEM - MEM BETWEEN FOODS
larv.ranef.wide.tax$diff.bmem.lmem <- larv.ranef.wide.tax$Bamboo_MEM - larv.ranef.wide.tax$Larval_MEM

##All & eco
##EMO - EMO BETWEEN FOODS
larv.ranef.wide.tax$diff.ball.lall <- larv.ranef.wide.tax$Bamboo_ALL - larv.ranef.wide.tax$Larval_ALL
##MEM - MEM BETWEEN FOODS
larv.ranef.wide.tax$diff.beco.leco <- larv.ranef.wide.tax$Bamboo_ECO - larv.ranef.wide.tax$Larval_ECO

##MMO - EMO BETWEEN FOODS
# larv.ranef.wide.tax$diff.bmmo.lemo <- larv.ranef.wide.tax$Larval_MMO - larv.ranef.wide.tax$Bamboo_EMO
# ##MMO - MEM BETWEEN FOODS
# larv.ranef.wide.tax$diff.bmmo.lmem <- larv.ranef.wide.tax$Larval_MMO - larv.ranef.wide.tax$Bamboo_MEM

##EMO - MMO BETWEEN FOODS
#larv.ranef.wide.tax$diff.lemo.bmmo <- larv.ranef.wide.tax$Larval_EMO - larv.ranef.wide.tax$Bamboo_MMO

##EMO - MEM BETWEEN FOODS
#larv.ranef.wide.tax$diff.lemo.bmem <- larv.ranef.wide.tax$Larval_EMO - larv.ranef.wide.tax$Bamboo_MEM

##MEM - MMO BETWEEN FOODS
#larv.ranef.wide.tax$diff.lmem.bmmo <- larv.ranef.wide.tax$Larval_MEM - larv.ranef.wide.tax$Bamboo_MMO
##MEM - EMO BETWEEN FOODS
#larv.ranef.wide.tax$diff.lmem.bemo <- larv.ranef.wide.tax$Larval_MEM - larv.ranef.wide.tax$Bamboo_EMO

3.1.2 Within foods

##MMO - EMO WITHIN BAMBOO
#larv.ranef.wide.tax$diff.bmmo.bemo <- larv.ranef.wide.tax$Bamboo_MMO - larv.ranef.wide.tax$Bamboo_EMO
##MMO - MEM WITHIN BAMBOO
larv.ranef.wide.tax$diff.bmmo.bmem <- larv.ranef.wide.tax$Bamboo_MMO - larv.ranef.wide.tax$Bamboo_MEM
##EMO - MEM WITHIN BAMBOO
larv.ranef.wide.tax$diff.bemo.bmem <- larv.ranef.wide.tax$Bamboo_EMO - larv.ranef.wide.tax$Bamboo_MEM

##MMO - EMO WITHIN LARVAL FOOD
#larv.ranef.wide.tax$diff.lmmo.lemo <- larv.ranef.wide.tax$Larval_MMO - larv.ranef.wide.tax$Larval_EMO
##MMO - MEM WITHIN LARVAL FOOD
larv.ranef.wide.tax$diff.lmmo.lmem <- larv.ranef.wide.tax$Larval_MMO - larv.ranef.wide.tax$Larval_MEM
##EMO - MEM WITHIN LARVAL FOOD
larv.ranef.wide.tax$diff.lemo.lmem <- larv.ranef.wide.tax$Larval_EMO - larv.ranef.wide.tax$Larval_MEM

3.2 Saving

#write.csv(larv.ranef.wide.tax,file="larv.ranef.wide.tax.all.csv")
#write.csv(larv.ranef.wide.tax,file="larv.ranef.wide.tax.csv")

3.3 Checking one out as an example

larv.ranef.diffs <- larv.ranef.wide.tax %>%
  select(id, starts_with("diff.")) %>%
  pivot_longer(
    cols = starts_with("diff."),
    names_to = "comparison",
    values_to = "difference"
  )

#larv.ranef.diffs.asv1 <- subset(larv.ranef.diffs,id=="ASV0001")

#ggplot(larv.ranef.diffs.asv1,aes(x=comparison,y=difference))+
#  geom_bar(stat="identity")+
#  theme(axis.text.x=element_text(angle=45,hjust=1))

##alternatively
ggplot(larv.ranef.wide.tax,aes(x=reorder(asv.gen,diff.bemo.lemo),y=diff.bemo.lemo))+
  geom_bar(stat="identity")+
  theme(axis.text.x=element_text(angle=45,hjust=1))+
  ggtitle("EMO - bamboo minus larval food")

##individual ASVs within treatments
ps.lar.trim.rel = transform_sample_counts(ps.lar.trim, function(x) x / sum(x))
ps.lar.trim.rel.trt <- subset_samples(ps.lar.trim.rel,Microbe_treatment=="MEM")
ps.lar.trim.rel.trt.asv <- subset_taxa(ps.lar.trim.rel.trt,id=="ASV0027")

plot_bar(ps.lar.trim.rel.trt.asv)+
  facet_wrap(Food_type~Microbe_treatment)

3.4 Plotting

larv.ranef.wide.tax <- read.csv("larv.ranef.wide.tax.csv")

3.4.1 Mosq. microbes

ggplot(larv.ranef.wide.tax,aes(x=reorder(asv.gen,diff.bmmo.lmmo),y=diff.bmmo.lmmo))+
  #geom_bar(stat="identity",position=position_dodge(1))+
  geom_bar(stat="identity")+
  #scale_fill_manual(values=c("#C6A947","#0085A6","#84D1B5"))+
  theme_cowplot()+
  theme(axis.text.x=element_text(angle=90,hjust=1))+
  ggtitle("MMO - bamboo minus larval food")

3.4.2 Env. microbes

ggplot(larv.ranef.wide.tax,aes(x=reorder(asv.gen,diff.bemo.lemo),y=diff.bemo.lemo))+
  #geom_bar(stat="identity",position=position_dodge(1))+
  geom_bar(stat="identity")+
  #scale_fill_manual(values=c("#C6A947","#0085A6","#84D1B5"))+
  theme_cowplot()+
  theme(axis.text.x=element_text(angle=90,hjust=1))+
  ggtitle("EMO - bamboo minus larval food")

3.4.3 Mosq + env microbes

ggplot(larv.ranef.wide.tax,aes(x=reorder(asv.gen,diff.bmem.lmem),y=diff.bmem.lmem))+
  #geom_bar(stat="identity",position=position_dodge(1))+
  geom_bar(stat="identity")+
  #scale_fill_manual(values=c("#C6A947","#0085A6","#84D1B5"))+
  theme_cowplot()+
  theme(axis.text.x=element_text(angle=90,hjust=1))+
  ggtitle("MEM - bamboo minus larval food")

3.5 Heat map

library("pheatmap")
library("RColorBrewer")

#diff.df <- larv.ranef.wide.tax[,c("asv.gen","diff.beco.leco","diff.bemo.lemo","diff.bmmo.lmmo","diff.bmem.lmem","diff.ball.lall")]
diff.df <- larv.ranef.wide.tax[,c("asv.gen","diff.bemo.lemo","diff.bmmo.lmmo","diff.bmem.lmem")]

row.names(diff.df) <- diff.df$asv.gen
diff.df2 <- diff.df[,2:4]
#diff.df2 <- diff.df[,2:6]
diff.mat <- as.matrix(diff.df2)

#row.names(otu.rel) == row.names(sam.rel)

pheatmap(diff.mat,
        color = colorRampPalette(c("#8C510A", "white", "darkolivegreen4"))(100),
        scale="row",
        cluster_rows=T,
        cluster_cols=T,
        #color = colorRampPalette(brewer.pal(n=7,name="BrBG"))(100)
        #color = colorRampPalette(brewer.pal(n=7,name="PiYG"))(100)
        #color = colorRampPalette(rev(brewer.pal(n=7,name="PiYG")))(100)
)

"#8C510A"
## [1] "#8C510A"

3.6 Plotting - all treatments

larv.ranef.wide.tax <- read.csv("larv.ranef.wide.tax.all.csv")

3.6.1 Mosq. microbes

ggplot(larv.ranef.wide.tax,aes(x=reorder(asv.gen,diff.bmmo.lmmo),y=diff.bmmo.lmmo))+
  #geom_bar(stat="identity",position=position_dodge(1))+
  geom_bar(stat="identity")+
  #scale_fill_manual(values=c("#C6A947","#0085A6","#84D1B5"))+
  theme_cowplot()+
  theme(axis.text.x=element_text(angle=90,hjust=1))+
  ggtitle("MMO - bamboo minus larval food")

3.6.2 Env. microbes

ggplot(larv.ranef.wide.tax,aes(x=reorder(asv.gen,diff.bemo.lemo),y=diff.bemo.lemo))+
  #geom_bar(stat="identity",position=position_dodge(1))+
  geom_bar(stat="identity")+
  #scale_fill_manual(values=c("#C6A947","#0085A6","#84D1B5"))+
  theme_cowplot()+
  theme(axis.text.x=element_text(angle=90,hjust=1))+
  ggtitle("EMO - bamboo minus larval food")

3.6.3 Mosq + env microbes

ggplot(larv.ranef.wide.tax,aes(x=reorder(asv.gen,diff.bmem.lmem),y=diff.bmem.lmem))+
  #geom_bar(stat="identity",position=position_dodge(1))+
  geom_bar(stat="identity")+
  #scale_fill_manual(values=c("#C6A947","#0085A6","#84D1B5"))+
  theme_cowplot()+
  theme(axis.text.x=element_text(angle=90,hjust=1))+
  ggtitle("MEM - bamboo minus larval food")

3.7 Heat map

library("pheatmap")
library("RColorBrewer")

diff.df <- larv.ranef.wide.tax[,c("asv.gen","diff.beco.leco","diff.bemo.lemo","diff.bmmo.lmmo","diff.bmem.lmem","diff.ball.lall")]
#diff.df <- larv.ranef.wide.tax[,c("asv.gen","diff.bemo.lemo","diff.bmmo.lmmo","diff.bmem.lmem")]

row.names(diff.df) <- diff.df$asv.gen
#diff.df2 <- diff.df[,2:4]
diff.df2 <- diff.df[,2:6]
diff.mat <- as.matrix(diff.df2)

#row.names(otu.rel) == row.names(sam.rel)

pheatmap(diff.mat,
        color = colorRampPalette(c("#8C510A", "white", "darkolivegreen4"))(100),
        scale="row",
        cluster_rows=T,
        cluster_cols=T,
        #color = colorRampPalette(brewer.pal(n=7,name="BrBG"))(100)
        #color = colorRampPalette(brewer.pal(n=7,name="PiYG"))(100)
        #color = colorRampPalette(rev(brewer.pal(n=7,name="PiYG")))(100)
)

"#8C510A"
## [1] "#8C510A"

4 Stopped updating script here

5 ggpredict

#install.packages("ggeffects")
library("ggeffects")

setwd("~/nicolagk@hawaii.edu - Google Drive/My Drive/Bamboo_mesos/Bamboo_mesos/04.mixed_models")

larv.mod.nb1 <- readRDS("larv.mod.nb1.rds")

#ps.trim.wtr <- readRDS("ps.exp.wtr.trim.rds")
#ps.trim.wtr
ps.trim.lar <- readRDS("ps.exp.lar.trim.rds")
ps.trim.lar
pred.lar.nb1 <- ggpredict(larv.mod.nb1, terms=c("variable", "Food_type","Microbe_treatment"), type="re")

plot(pred.lar.nb1)+
  theme(axis.text.x=element_text(angle=90))+
  theme_classic()

# pred.lar.nb2 <- ggpredict(larv.mod.nb2, terms=c("variable", "Food_type","Microbe_treatment"), type="re")
# plot(pred.lar.nb2)

tax.trim.lar <- data.frame(ps.trim.lar@tax_table)

##make same column names
pred.lar.nb1$id <- pred.lar.nb1$x
pred.lar.tax <- merge(tax.trim.lar,pred.lar.nb1,by="id")

ggplot(pred.lar.tax,aes(x=group,y=predicted,color=facet))+
  geom_point()+
  facet_wrap(~facet)

pred.lar.tax.mmo <- subset(pred.lar.tax,facet=="MMO")

ggplot(pred.lar.tax.mmo,aes(x=group,y=predicted,fill=id))+
  #geom_bar(stat="identity",position=position_dodge(0.5))+
  geom_bar(stat="identity")+
  facet_wrap(~Genus,scales="free")+
  theme(legend.position="none")

pred.lar.tax.bamb <- subset(pred.lar.tax,group=="Bamboo")

ggplot(pred.lar.tax.bamb,aes(x=facet,y=predicted,fill=id))+
  geom_bar(stat="identity")+
  facet_wrap(~Family,scales="free")+
  theme(legend.position="none")

Trying some different plots

library("cowplot")

pred.lar.tax.comp <- subset(pred.lar.tax,facet!="ECO"&facet!="ALL")

pred.lar.tax.comp.bam <- subset(pred.lar.tax.comp,group=="Bamboo")

ggplot(pred.lar.tax.comp.bam,aes(x=facet,y=predicted,fill=facet))+
  facet_wrap(~id,scales="free")+
  geom_bar(stat="identity")+
  theme_cowplot()

##heatmap
library("pheatmap")
library("RColorBrewer")

library("tidyr")

#pred.lar.tax.comp$asv_food <- paste0(pred.lar.tax.comp$id,"_",pred.lar.tax.comp$group)

colnames(pred.lar.tax.comp)

##remove not needed columns
pred.lar.min <- pred.lar.tax.comp %>%
  select(predicted,facet,group,id)

##bamboo
pred.lar.min.bam <- subset(pred.lar.min,group=="Bamboo")

pred.lar.min.bamw <- pred.lar.min.bam %>%
  pivot_wider(names_from=id,values_from=predicted) %>%
  as.data.frame()

row.names(pred.lar.min.bamw) <- pred.lar.min.bamw$facet

pred.lar.min.bamw2 <- pred.lar.min.bamw[,3:ncol(pred.lar.min.bamw)]

pheatmap(as.matrix(pred.lar.min.bamw2),
        color = colorRampPalette(c("white", "orange", "red"))(100),
        #color = colorRampPalette(rev(brewer.pal(n=7,name="PiYG")))(100),
        #scale="column",
        cluster_cols = F,
        cluster_rows=F
)

##larval food
pred.lar.min.lar <- subset(pred.lar.min,group=="Larval food")

pred.lar.min.larw <- pred.lar.min.lar %>%
  pivot_wider(names_from=id,values_from=predicted) %>%
  as.data.frame()

row.names(pred.lar.min.larw) <- pred.lar.min.larw$facet

pred.lar.min.larw2 <- pred.lar.min.larw[,3:ncol(pred.lar.min.larw)]

pheatmap(as.matrix(pred.lar.min.larw2),
        color = colorRampPalette(c("white", "orange", "red"))(100),
        #color = colorRampPalette(rev(brewer.pal(n=7,name="PiYG")))(100),
        #scale="column",
        cluster_cols = F,
        cluster_rows=F
)

library("ggrepel")

##correlating what happened in bamboo & larval food
pred.lar.min$group <- gsub("Larval food","Larval",pred.lar.min$group)

pred.lar.min.cor <- pred.lar.min %>%
  pivot_wider(names_from=group,values_from=predicted)

pred.lar.min.cor$facet <- factor(pred.lar.min.cor$facet,levels=c("MMO","EMO","MEM"))

##full plot
ggplot(pred.lar.min.cor,aes(x=Bamboo,y=Larval,color=facet,shape=facet,label=id))+
  geom_point(size=2)+
  geom_abline(intercept=0,slope=1,linetype="dotted")+
  geom_text_repel(color="black")+
  theme_cowplot()+
  ylab("Standard diet")+
  xlab("Bamboo leaves")+
  scale_color_manual(values=c("#C6A947","#84D1B5","#0085A6"))

##zooming into the low section
ggplot(pred.lar.min.cor,aes(x=Bamboo,y=Larval,color=facet,shape=facet,label=id))+
  geom_point()+
  geom_abline(intercept=0,slope=1,linetype="dotted")+
  geom_abline(intercept=0,slope=0.5,linetype="dotted")+
  geom_abline(intercept=0,slope=2,linetype="dotted")+
  xlim(0,100)+
  ylim(0,100)+
  geom_text_repel(color="black")+
  theme_cowplot()+
  ylab("Standard diet")+
  xlab("Bamboo leaves")+
  scale_color_manual(values=c("#C6A947","#84D1B5","#0085A6"))

mean(sample_sums(ps.lar.trim))
#47762.4
1000/47762.4

##max 2k - larval food
ggplot(pred.lar.min.cor,aes(x=Bamboo,y=Larval,color=facet,shape=facet,label=id))+
  geom_point()+
  geom_abline(intercept=0,slope=1,linetype="dotted")+
  geom_abline(intercept=0,slope=0.5,linetype="dotted")+
  geom_abline(intercept=0,slope=2,linetype="dotted")+
  #xlim(0,4000)+
  ylim(1000,30000)+
  geom_text_repel(color="black")+
  theme_cowplot()+
  ylab("Standard diet")+
  xlab("Bamboo leaves")+
  scale_color_manual(values=c("#C6A947","#84D1B5","#0085A6"))

##max 2k - bamboo
ggplot(pred.lar.min.cor,aes(x=Bamboo,y=Larval,color=facet,shape=facet,label=id))+
  geom_point()+
  geom_abline(intercept=0,slope=1,linetype="dotted")+
  geom_abline(intercept=0,slope=0.5,linetype="dotted")+
  geom_abline(intercept=0,slope=2,linetype="dotted")+
  xlim(1000,30000)+
  #ylim(1000,30000)+
  geom_text_repel(color="black")+
  theme_cowplot()+
  ylab("Standard diet")+
  xlab("Bamboo leaves")+
  scale_color_manual(values=c("#C6A947","#84D1B5","#0085A6"))

Checking what the actual data look like

ps.trim.lar.comp <- subset_samples(ps.lar.trim,Microbe_treatment!="ALL"&Microbe_treatment!="ECO")
ps.trim.lar.comp

ps.lar.trim.rel = transform_sample_counts(ps.trim.lar.comp, function(x) x / sum(x))
ps.lar.trim.rel

ps.asv1 <- subset_taxa(ps.trim.lar.comp,id=="ASV0001")
ps.asv1@otu_table
plot_bar(ps.asv1,x="Sample")+
  facet_wrap(Microbe_treatment~Food_type)
plot_bar(ps.asv1)+
  facet_wrap(Microbe_treatment~Food_type,scales="free")

ps.asv3 <- subset_taxa(ps.lar.trim.rel,id=="ASV0003")
ps.asv3
plot_bar(ps.asv3,x="Microbe_treatment")+
  facet_wrap(~Food_type)
plot_bar(ps.asv3)+
  facet_wrap(Microbe_treatment~Food_type,scales="free")

ps.asv4 <- subset_taxa(ps.lar.trim.rel,id=="ASV0004")
ps.asv4
plot_bar(ps.asv4,x="Microbe_treatment")+
  facet_wrap(~Food_type)

ps.asv6 <- subset_taxa(ps.lar.trim.rel,id=="ASV0006")
ps.asv6
plot_bar(ps.asv6,x="Microbe_treatment")+
  facet_wrap(~Food_type)

ps.asv7 <- subset_taxa(ps.lar.trim.rel,id=="ASV0007")
ps.asv7
plot_bar(ps.asv7,x="Microbe_treatment")+
  facet_wrap(~Food_type)

ps.asv19 <- subset_taxa(ps.trim.lar.comp,id=="ASV0019")
ps.asv19
plot_bar(ps.asv19,x="Microbe_treatment")+
  facet_wrap(~Food_type)
plot_bar(ps.asv19)+
  facet_wrap(Microbe_treatment~Food_type,scales="free")

ps.asv33 <- subset_taxa(ps.trim.lar.comp,id=="ASV0033")
ps.asv33
plot_bar(ps.asv33,x="Microbe_treatment")+
  facet_wrap(~Food_type)
plot_bar(ps.asv33)+
  facet_wrap(Microbe_treatment~Food_type,scales="free")

ps.asv54 <- subset_taxa(ps.lar.trim.rel,id=="ASV0054")
ps.asv54
plot_bar(ps.asv54,x="Microbe_treatment")+
  facet_wrap(~Food_type)
plot_bar(ps.asv54)+
  facet_wrap(Microbe_treatment~Food_type)

larv.ranef.wide.tax[6,]

         
#         )+
#  facet_wrap(Food_type~Microbe_treatment)